Travelling waves and impact-parameter correlations 
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It is usually assumed that the high-energy evolution of partons in QCD remains local in coordinate 
space. In particular, fixed impact-parameter scattering is thought to be in the universality class of 
one-dimensional reaction-diff'usion processes as if the evolutions at different points in the transverse 
plane became uncorrelated through rapidity evolution. We check this assumption by numerically 
comparing a toy model with QCD-like impact-parameter dependence to its exact counterpart with 
uniform evolution in impact-parameter space. We find quantitative differences, but which seem to 
amount to a mere rescaling of the strong coupling constant. Since the rescaling factor does not 
show any strong Qs-dependence, we conclude that locality is well verified, up to subleading terms 
at small as- 



I. INTRODUCTION 

Recently, equations describing evolution towards high 
energy in QCD, including saturation effects [iIjSSIJjS 
Mi , have been proven [7] to belong to the universality class 
of the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) 
equation. It has also been realised that additional 
contributions were to be included in order to fully satisfy 
unitarity, in particular, in a way consistent with boost- 
invariance. Though only some partial results exist for 
a more complete set of evolution equations [^, [l3, fl]| . 
the evolution of partonic states and their scattering at 
high energy is conjectured [l2| to be a process of the 
reaction-diffusion type. The QCD evolution equations 
are thus thought to belong to the universality class of 
the stochastic FKPP (sFKPP) equation [13] which reads 

dtu{t, x) — d'^u{t, x) + u(t, x) — u^{t, x) 



+ ey/2{u{t, x) - u'^(t, x)) v{t, x) (1) 

where v \s a normal Gaussian noise uncorrelated in space 
and time; t corresponds to the rapidity variable in QCD 
(or the logarithm of the energy), and x is the logarithm of 
the transverse scale (momentum or distance) relevant to 
characterise the partons in the considered process, while 
e, the strength of the noise, is of the order of the strong 
coupling constant a^; u can be thought of an event-by- 
event scattering amplitude, whose average over events is 
the physical amplitude. Such an equation is supposed to 
hold at a specific impact parameter. 

A convenient picture for high-energy scattering is the 
colour dipole model [ij], whose numerical implemen- 
tation in terms of Monte-Carlo code event generators 
[isl Ha . \Vn . [l8j has been particularly useful to investi- 
gate saturation issues in QCD. The picture is the fol- 
lowing: a hadron may be represented as a set of dipoles 
at the time of the interaction, which is constructed by 
successive dipole splittings (i.e. emissions of giuons in 
the large- A^c approximation) and mergings (or any other 
nonlinear process that would tame the growth of their 
number) whose probabilities are computed from QCD 



vertices. The set of dipoles resulting from this process 
eventually interacts with a target by exchanging giuons. 

One condition for the one-dimensional^ sFKPP equa- 
tion to be relevant to high energy scattering is that the 
evolution remains local in impact parameter: indeed, if 
this was not the case, then correlations between different 
points in impact-parameter space would invalidate one- 
dimensional equations like ([!]) in QCD. At first sight, it 
is not at all obvious that this condition should hold: the 
dipole evolution kernel is singular and allows for split- 
tings into areas that lie very far from each other in co- 
ordinate space. This nonlocal behaviour may a priori 
cause migrations of the dipole chains of successive split- 
tings over large distances in the transverse plane, which 
in turn may induce correlations between different impact 
parameters. 

In Ref. [l^l, this independence of the different im- 
pact parameters was assumed on the basis of a simple 
analytical estimate of the mean number of dipoles pro- 
duced coUinearly. Capturing the full complexity of the 
QCD evolution to achieve a more precise understanding 
still seems out of reach of an analytical approach. On 
the other hand, to our knowledge, all recent numerical 
studies of high-energy scattering at saturation (including 
gluon-number fluctuations) and its relation to the sFKPP 
equation (see e.g. [l3,[20|) relied on the one-dimensional 
approximation. 

Our goal in this paper is to test the validity of a 
one-dimensional formulation such as Eq. ([1]) in a model 
that has impact-parameter dependence and that we shall 
study numerically. We could consider the full colour 
dipole model and use the already-existing Monte-Carlo 
code in Ref. [l6j , but its non trivial 2-dimensional impact- 
parameter structure makes it far too complex to handle, 
especially since we want to run the evolution over a time 



^ Through this paper, by one- dimensional we refer to models with 
a single spatial dimension {e.g. the dipole size), in addition to 
the evolution variable {e.g. time or rapidity). 



long enough to observe the asymptotic properties. In- 
stead, we shall build a simpler model, with discretised 
dipole sizes and a single dimension of impact parameter, 
which reproduces the main features of full QCD as far as 
the singularities are concerned. 

The structure of the paper is as follows. Sec. |TT] re- 
calls the essential features of the colour dipole model 
and its interpretation in a statistical mechanics language. 
In Sec. mil we construct a toy model which incorporates 
an impact-parameter dependence, at variance with ex- 
isting models of QCD evolution. We also introduce a 
corresponding model without impact parameter which 
would be fully equivalent to the first one in the absence 
of correlations in impact parameter. Finally, in Sec. IIVI 
we present the numerical results on travelling waves and 
their impact-parameter correlations obtained within this 
model and compare with the model without impact pa- 
rameter dependence. 



II. ANALYTICAL PRELIMINARIES 

A. Scattering in the dipole model 

In the colour dipole model [lj|, high-energy evolution 
is viewed in the following way. Two hadrons (which are 
asymptotically sets of colourless qq dipoles, for the sim- 
plicity of the argument) develop under rapidity evolu- 
tion highly occupied Fock-states which themselves may 
be seen as collections of colour dipoles. The building 
up of the states of each hadron is specified by providing 
the splitting rate of a dipole whose endpoints have trans- 
verse coordinates (a;o,a;i) into two dipoles (2;o,a;2) and 
(xi,a:2) as the result of a gluon emission at position X2- 
It reads 111 



dP 

d{ay) 



(xoi -^ a;o2,a;i2) 



X, 



01 



X2 



^02^^12 



27r 



(2) 



This splitting process is actually supplemented by a sat- 
uration mechanism, which limits the local density of 
dipoles to about l/ofj. Unfortunately, it has not been 
formulated in the dipole model. (For a study of the prob- 
lems that one has to face to incorporate saturation effects 
in this framework, see e.g. Ref. [2l|]). It is sometimes 
seen as a progressive slowing down of the dipole emis- 
sion rate once the local density is of the order of l/a^ 
[3, [22|. It could also be recombination of dipoles, which 
corresponds to the reabsorption of gluons [2l|, or some 
form of colour reconnection ("swings") among the con- 
stituent gluons of the dipoles [iM| . But the very nature 
of the mechanism is not crucial to our discussion, and 
when building our toy model, we pick the simplest one 
to implement numerically. 

At the time of the interaction, pairs of dipoles made of 
one dipole from the left-moving hadron and one dipole 
from the right- moving hadron, may exchange gluons, pro- 
vided that they are of similar sizes and sit at the same 



impact parameter. It is actually simpler to think of the 
process in the rest frame of one of the hadrons: then, 
only the moving hadron (the probe) undergoes evolu- 
tion, while the other one (the target) remains in its bare 
asymptotic state. 

We need to clarify what is meant by two dipoles sitting 
at the same impact parameter (up to a distance of the 
order of their size). When a highly evolved hadron is 
probed by an elementary dipole of size xqi — xq — xi, 
the scattering amplitude T is related to the density n in 
the following way: 



T{y,xoi,2^^)^ 

2 ko 



27r 277 



X los 



zi?\xi - za\ 



\xo -zopki - zi\- 



■n(y,^oi,i4^), (3) 



where the summation goes over the positions of the 
dipoles in the Fock state of the evolved hadron. The 
first argument of T is the rapidity, the second one the 
dipole size and the last one the impact parameter. This 
formula means that T roughly counts the dipoles which 
have a size of the order of |xoi| and sit in a region of 
impact parameter delimited by a circle of radius |xoi| 
around the position (or impact parameter) [xq + a;i)/2. 
From the dipole splitting probability ([2]), we can infer 
the evolution equation for T . The simplest version of this 
equation is obtained when the target is a very large nu- 
cleus, that is, made of a large (formally infinite) number 
of bare dipoles. One gets 



aa„r(y,xoi,^4^) 



d^X2 Xqi 
ZTT a^02-''12 L 



Tiy, xo2, 



gQ+^2 ' 



+ T(y, X12, ^^4^) - T{y, xqi, ^^4^) 

- T{y, a;o2, ^^4^)r(y, a:i2, ^^4^) 



(4) 



This is the well-known Balitsky-Kovchegov (BK) equa- 
tion [3, 0] . Its linear part is the Balitsky-Fadin-Kuraev- 
Lipatov (BFKL) equation [23|, while its nonlinear part 
accounts for multiple scattering on the dense target. For 
a general target, additional unitarity effects must be 
taken into account [Sj, [12|] , and the corresponding equa- 
tion is not yet fully known, though some partial results 
covering the main features are available [3, [l^, \lM, [2J] ■ 
For example, at large- iVc and in the two-gluon-exchange 
approximation, it has been shown [^ that ([4]) has to be 
supplemented by a non-local noise term taking into ac- 
count the gluon-number fluctuations when the target is 
dilute. 



B. Travelling waves 

Recent research has established a link between high 
energy QCD evolution and reaction-diffusion processes. 
Indeed, the Balitsky-Kovchegov equation at fixed impact 



parameter, in the diffusive approximation^, was sfiown 
to be identical to the FKPP equation \^. Moreover, 
as recalled in the Introduction, it is thought that full 
high-energy QCD beyond the approximations assumed 
to establish the BK equation lies in the universality class 
of reaction-diffusion processes [i2[. The stochastic ver- 
sion of the FKPP equation (see Eq. ([1])) then becomes 
relevant, with a noise term of the order of the strong 
coupling constant as- The interesting outcome of this 
correspondence consists in universal results for the scat- 
tering amplitudes that can be analytically computed in 
the asymptotic regime of very small as and large rapidi- 
ties though, for larger values of the couphng, the same 
behaviours are observed numerically [iSl, [20, l2^ ■ Even if 
the precise form of the evolution equations are not fully 
known, all the models presented so far have shown the 
same asymptotic properties. 

Generally speaking, equations of the form ([4]) admit 
travelling wave solutions in the variable log(l/a::§]^) when 
the dependency upon the impact parameter is ignored. 
The main features of these waves only depend on the 
Mellin moments xil) of the splitting probability ([2]). 
These solutions are attractors in the sense that a large 
class of initial conditions (actually all initial conditions 
that are physical for the considered processes) converge 
to them at large rapidity. The position Ps{y) = log Q^{y) 
of the wave defines the saturation scale Qs{y)- Thus, 
the saturation scale and its rapidity dependence emerge 
naturally from the travelling-wave solutions. At large 
rapidities, the speed of the wave dayPs{y) is [1, 25,, 26|] 



where 7c satisfies 
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(5) 



(6) 



The wave front decays exponentially in the sinall-|xoi| 
region like 



T ^ exp -7c lo; 



1 



m 



(7) 



Geometric scaling [27| is the physical phenomenon that 
corresponds to the fact that the amplitude only depends 
on the combination XQiQ^{y) rather than on each vari- 
able separately. 

If stochasticity is taken into account, then the proper- 
ties of the solutions are modified. The travelling waves 
become noisy, corresponding to the fact that each indi- 
vidual event is made of a discrete number of dipoles in- 
troducing fluctuations important in the tail of the front 



where the dipole density is of order 1 (or, T ^ a^). The 
logarithm of the squared saturation scale acquires a vari- 
ance that grows linearly with rapidity. The velocity of 
the fronts, that is to say, the rates of growth of log Q^{y) 
averaged over the realisations, reach an asymptotic value 
smaller than that coming out of the BK equation by a 
shift proportional to 1/ log^(l/a^). The physical scatter- 
ing amplitude is obtained by averaging over events, and 
gets tilted [2^ with respect to what it would be if the 
BK equation held. Geometric scaling is broken \§\ and 
replaced by diffusive scaling @, M, \1M '■ 



(T) = erfc 



^ay/log\l/al) 



(8) 



^ The diffusive approximation corresponds to expanding the BFKL 
kernel to second order. 



which reflects the fact that the variance of the front po- 
sition scales like ay / \o^ {1 / aj.) . 

It is important to note that the wave velocity gets 
closer to the mean-field result ([5|) when as decreases, 
and not surprisingly, the variance of the front position 
gets smaller in the same limit, which is the sign that 
fluctuations get milder. This limit corresponds to large 
maximum local densities of dipoles, for which statistical 
fluctuations are indeed expected to vanish, and thus a 
mean-field approximation may indeed be justified. 

These results are believed to be rigorously true for one- 
dimensional models of reaction-diffusion type. Their ex- 
tension to QCD actually relies in particular on the postu- 
late that the impact-parameter variable does not play a 
role in equations such as dU and in their stochastic exten- 
sions, for the latter should be "local" enough. However, 
this fundamental assumption has never been checked, 
and this is precisely what we intend to do numerically 
in this paper. Before we turn to this task (Sec. IIII[) . let 
us recall the arguments in favour of the decoupling of the 
travelling waves at different impact parameters. 



C. Picture of the evolution including impact 
parameter 

Let us start with a single dipole at rest, and bring 
it gradually to a higher rapidity. During this process, 
this dipole may be replaced by two new dipoles, which 
themselves may split, and so on, eventually producing a 
chain of dipoles. Figure [1] pictures one realization of such 
a chain. 

According to Eq. ([2]), splittings to smaller-size dipoles 
are favoured, and thus, one expects that the sizes of the 
dipoles get smaller on the average, and that in turn, 
the successive splittings become more local. The dipoles 
around region "1" and those around region "2" should 
have an independent evolution beyond the stage pictured 
in the figure: further splittings will not mix in impact 
parameter space, and thus, the travelling waves around 
these regions should be uncorrelated. For a dipole in re- 
gion 1 of size r to migrate to region 2, it should first 




FIG. 1: Sample dipole configuration in impact-parameter 
space. The points represent the gluons, and the fine is the 
dipole chain after evolution. The initial dipole is represented 
by a dashed line, and the first splitting by two dotted lines. 
For the two impact parameters denoted by a "+" sign (la- 
belled 1 and 2) their further respective evolutions are be- 
lieved to be uncorrelated. The dipoles whose both endpoints 
are empty circles are those which are seen at these peculiar 
impact parameters. For clarity, we have omitted the large 
number of tiny dipoles produced at every endpoints. 



split into a dipole whose size is of the order of the dis- 
tance A6 between regions 1 and 2, up to a multiplicative 
uncertainty of order 1. (We assume in this discussion 
that the dipoles in region 2 relevant to the propagation 
of the local travelling waves, that is, those which are in 
the bulk of the wave front, also have sizes of order r). 
Roughly speaking, the rate of such splittings may be es- 
timated from the dipole splitting probability ([2]): it is 
of order a(r^/(A6)^)^, while the rate of splittings of the 
same dipole into a dipole of similar size in region 1 is of 
order a. Thus the first process is strongly suppressed as 
soon as regions 1 and 2 are more distant than a few units 
of r. Note that for A6 > l/Qs, saturation may further 
reduce the emission of the first, large, dipole leading to 
an even stronger suppression of the estimated rate. 

There is a second case that we should worry about. 
What could also happen is that some larger dipole has, 
by chance, one of its endpoints tuned to the vicinity of 
the coordinate one is looking at (at a distance which is 
at most |Ar| ^ l/Qs{Y)), and easily produces a large 
number of dipoles there. In this case, the position of the 
travelling wave at that impact parameter would suddenly 
jump. If such events were frequent enough, then they 
would modify the average wave velocity and thus the one- 
dimensional picture. We may give a rough estimate of 
the rate at which dipoles of size smaller than Ar are 
produced. Assuming local uniformity for the distribution 
of the emitting dipoles, the rate (per unit of ay) of such 
events can be written 



rfVo 



ro>Ar 



d en(ro) 

e<Ar Vo 



1 



27r £2(ro — e)' 



:, (9) 



where we integrate over large dipoles of size tq > Ar 
emitting smaller dipoles (of size e < Ar) with a proba- 
bihty d^ero/(27re^(ro— e)^). The factor (e/rg)^ accounts 
for the fact that one endpoint of the dipole of size ro has 



to be in a given region of size e in order to emit the 
dipoles at the right impact parameter. To estimate this 
expression, we first use n{ro) — T{ro)/a'^ and use for T 
the simplified expression 

T{ro) - e(ro - l/Qs) + {rlQly- e(l/Q, - ro), 

which splits the front into a saturated region (ro > l/Qs) 
and a tail with geometric scaling (ro < l/Qs)- Using 
ro — e ~ ro in the emission kernel, the integration is then 
easily performed and one finds a rate whose dominant 
term is 



2\7c 



TT iiAr)'Qi) 



2a2 



1 



■7c 



(10) 



For (Ar)2 < (al)^^'^" /Ql, that is, ahead of the bulk of 
the front, this term is parametrically less than 1 and is 
in fact of the order of the probability to find an object 
in this region that contributes to the normal evolution 
of the front [23]. Hence there is no extra contribution 
due to the fact that there are many dipoles around at 
different impact parameters. 

However, it seems difficult to perform a precise cal- 
culation of these effects beyond order-of-magnitude esti- 
mates, and anyway, this average picture may be spoiled 
by the statistical event-by-event fluctuations. So, we do 
not have any insight on how to estimate more precisely 
this effect analytically. It is thus important to perform a 
numerical simulation and check the correlations of trav- 
elling waves at different impact parameters. 



III. CONSTRUCTION OF THE MODEL 

A. Simplified model with impact-parameter 
dependence 

1. Fock state evolution 

In order to arrive at a model that is tractable numer- 
ically, we only keep one transverse dimension instead of 
two in 3-|-l-dimensional QCD. However, we cannot con- 
sider genuine 2-|-l-dimensional QCD because we do not 
wish to give up the logarithmic coUinear singularities at 
X2 — xq and X2 = xi. Moreover, QCD with one dimen- 
sion less has very different properties at high energies 
[30[. A splitting rate which complies with our require- 
ments is: 



1 



Foil 
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d{ay) 4|a;o2||a;i2 



■dx-2 



(11) 



We can further simplify this probability distribution by 
keeping only its coUinear and infrared asymptotics (as in 
[3l|). If \xq2\ ^ |a^oi| (or the symmetrical case |xi2| ^ 
jsoil), the probability reduces to dx2/\x{)2\ {dx2/\xi2\ 
resp.). The result of the splitting is a small dipole {xq,X2) 
together with one close in size to the parent. So for sim- 
plicity we will just add the small dipole to the system 



and leave the parent unchanged. In the infrared region, 
a dipole of size |a;o2| 2> |a;oi| is emitted with a rate given 
by the large- |a;o2| limit of the above probability. The 
probability laws (I1]),((TT|) imply that a second dipole of 
similar size should be produced while the parent dipole 
disappears. To retain a behaviour as close as possible to 
that in the coUinear limit, we will instead just generate 
a single large dipole and maintain the parent. To do this 
consistently one must include a factor of two in the in- 
frared splitting rate, so as not to modify the average rate 
of production of large dipoles. 

In formulating our model precisely, let us focus first on 
the distribution of the sizes of the participating dipoles. 
(The simplifying assumptions made above enable one to 
choose the sizes and the impact parameters of the dipoles 
successively). We call r the modulus of the emitted 
dipole. To the modulus of its parent and Y — ay. The 
splitting rate pT|) reads, in this simplified model 



dPr, 



dY 



= y{r - roj 



r^dr 



0{ro 



dr 



(12) 



and the original parent dipole is kept. Logarithmic vari- 
ables are the relevant ones here, so we introduce 



p = logB(l/r) 



r = B- 



(13) 



where the base B will later be set to two in the actual 
numerical calculations. We can thus rewrite the dipole 
creation rate as 



dP, 



""^^ = e{po - p) BP-P° log B dp 

+ eip- po) log B dp. (14) 



dY 



To further simplify the model, we discretise the dipole 
sizes and consider a lattice in p with lattice spacing A 
(which will later be set to one) . This amounts to restrict- 
ing the dipole sizes to negative integer powers of B . The 
probability that a dipole at lattice site i {i.e. a dipole of 
size B^^^) creates a new dipole at lattice site j is 



dP,. 



r«+i dP, 



dY 




-i)A 



j > i 
j < i 



(15) 
(16) 



The rates dPi±/dY for a dipole at lattice site i to split 
to any lattice site j > i or j < i respectively are then 
given by 



dY 
dY 



dY 



L-1 

E 

'"^ dP 



A\ogB{L-i), (17) 



j=o 



dY 



= 1- B 



-iA 



(18) 



where we have restricted the lattice to < i < L, for 
obvious reasons related to the numerical implementation. 



Now we have to address the question of the impact 
parameter of the emitted dipole. In QCD, the coUinear 
dipoles are produced near the endpoints of the parent 
dipoles. Let us take a parent of size rp at impact param- 
eter bo. We set the emitted dipole (size r) at the impact 
parameter b such that 



6o± 



vq ±r X s 



(19) 



where s has uniform probability between and 1. It is 
introduced to obtain a continuous distribution of the im- 
pact parameter unaffected by the discretisation of r. This 
prescription is quite arbitrary in its details, but the lat- 
ter do not influence significantly the physical observables. 
Each of the two signs that appear in the above expres- 
sion is chosen to be either -I- or — with equal weights. 
We apply the same prescription when the emitted dipole 
is larger than its parent. 



2. Scattering amplitude 

We explained before that in QCD, the scattering am- 
plitude of an elementary probe dipole of size r^ = B^^^ 
with a dipole in an evolved Fock state is proportional to 
the number of objects which have a size of the same order 
of magnitude and which sit in an area of radius of order 
ri around the impact point of the probe dipole. Since 
in our case, the sizes are discrete, the amplitude is just 
given, up to a factor, by the number of dipoles that are 
exactly in the same bin of size as the probe, namely 

T{i, bo) = (a^/A) X #{dipoles of size B"^^ 
at impact parameter 6 such that \b — bo\ < fi/2}. (20) 



3. Saturation 

We now have to enforce unitarity, that is the condition 

T{t,b)<l (21) 

for any i and b. This condition is expected to hold due 
to gluon saturation in QCD. However, saturation is not 
included in the original dipole model. Nevertheless, as 
argued in Section |TT1 the asymptotic properties are not 
affected by the details on how we implement the condi- 
tion ([2T|) . The simplest choice is to veto splittings that 
would locally drive the amplitude to values larger than 
1. In practice, for each splitting that gives birth to a new 
dipole of size i at impact parameter 6, we compute T{i, b) 
and T(i,bi: ?'i/2), and throw away the produced dipole 
whenever one of these numbers gets larger than one. 

Given the definition of the amplitude T, this saturation 
rule implies that there is a maximum number of objects 
in each bin of size and at each impact parameter, which 
is equal to iVgat = A/a^. 



4- Note on the implementation 

The model is now completely specified by 
Eqs. (USD, dUl), dSOD and dUD- 

The implementation of the dipole splittings is quite 
straightforward, since the distribution of the sizes and 
impact parameters of the produced dipoles are very sim- 
ple. At a given point along the rapidity evolution, we 
choose the dipole that is going to split next. The dipole 
lifetimes only depend on their sizes, and are simply given 
by the inverse rates d{Pi+ + Pi-)/dY in Eq. ^TE\\ . Once 
this dipole has been chosen, we determine its size and 
impact parameter according to Eqs. pB)) . p^ . Next, we 
check that the conditions for the amplitude not to violate 
unitarity are satisfied. If this is the case, the new dipole 
is integrated in the Fock state. Otherwise, the dipole is 
rejected. 

The structure of the data has to be chosen carefully to 
save CPU time. Indeed, in order to compute the cross 
section T, one needs to search for all dipoles of a given 
size within a specified range of impact parameter about 
the probe dipole. Given that the number of dipoles grows 
very quickly with rapidity (a priori exponentially), it is 
crucial to be able to search these dipoles in a time that 
does not grow faster than the logarithm of their num- 
ber. To meet this requirement, an array indexed by the 
discrete logarithmic size i of the dipoles points to binary 
trees of dipoles of identical sizes, ordered by their impact 
parameter. 

The number of objects grows very fast due to the large 
number of tiny dipoles that may easily be produced, 
and consequently, the computer memory would be sat- 
urated quite early in the evolution. Imposing the above- 
mentioned cutoff L on the logarithmic size is not efficient 
enough to limit the number of objects. The reason is that 
L has to be taken quite high (we chose 50): if it is lower 
than that, the evolution gets hampered within the range 
of rapidity that is of interest for investigating saturation. 
Hence we decided to enforce a lower cutoff also on the ra- 
tio K = B~'^/\b — BnI of the size of each produced dipole 
to the distance between its impact parameter b and any 
impact parameter bN at which we would measure cross 
sections. Imposing this cutoff, we anticipate on the fact 
that a small dipole very far from the impact parameters 
of interest has only a tiny probability to migrate back 
through its further splittings; k has of course to be taken 
sufficiently small, and varied in order to evaluate its ef- 
fect. 

Restrictions due to machine accuracy must also be ad- 
dressed. Indeed, when dipoles become small through evo- 
lution, one must be able to resolve numerically equally 
small distances in impact parameter. Since absolute im- 
pact parameters will be of the order of 1 and are coded 
on 53 significant bits (typical double precision floating 
point numbers) , we have to limit the dipole sizes to about 
i < 53 log S/ log 2. A way to overcome this limitation is 
to use arbitrary precision arithmetic libraries |32j], but 



this is both more memory and CPU costly so we tried to 
avoid this solution. 



B. Fixed impact-parameter version 

We may now consider a similar model, but in which 
there is no impact-parameter dependence. We will call it 
"FIP" for Fixed Impact Parameter, while the complete 
model will be termed "AIP" (Any Impact Parameter). 



1. Formulation 

One may understand how to "remove" the impact pa- 
rameter from the AIP model as follows. The essential 
point is that the interaction between the probe and tar- 
get is local in impact parameter. That means that if a 
probe of size r is to interact with a dipole of size r in the 
target, then the probe should also be roughly within a 
transverse distance r of the target dipole. When explicit 
impact parameter information is maintained, one simply 
enforces this in one's determination of probe-target in- 
teractions. If instead one discards the impact parameter 
information, then one must find another way of account- 
ing for the fact that most target dipoles of size r will be 
too far from the probe to have a significant interaction 
with it. The simplest is to observe that the probability of 
a target dipole being close enough to the probe is r/ro, 
where rg is the initial size of the system. Therefore the 
number of dipoles of size r sufficiently close to a specific 
fixed impact parameter, n^-^'{r), is given in terms of the 
total number of dipoles of size r, n(r), by the relation 
rv-f'{r) — r/ron{r). One may reinterpret this in terms 
of an FlP-specific branching probability, which relates to 
the normal branching probability via 



dP. 



-^3 



fj dPi^j 



dY 



n dY 
1 (B^ 1) 



^(i-7)A dPj^j 



dY 

3 > i 
3 < i 



(22) 
(23) 



It is straightforward to verify that this reproduces the 
condition 'n!'^\r) = r/rori{r), and from eqs. (22,23) one 
then derives 



dP, 



if) 



i+ 



L-1 ^p(f) 



dY 



E 



1 _ R-(i-i)A 

'^' =AlogB^ ^^^^, (24) 



dY 
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dY 



E' 
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dY 



{B'^ -l)i. 



(25) 



Of course, the saturation rule in the FIP model is that 
a dipole is not produced if the number of objects of the 
same size is already larger than A'sat- 



S. Numerical implementation 

The FIP model is quite straightforward to implement. 
The splitting probabilities are as simple as those in the 
AIP model. Since again the dipoles have discrete sizes, it 
is enough to update an array indexed by the logarithmic 
dipolc size i, each cell of which contains the correspond- 
ing number of dipoles. 



the value of A. From the expression Eq. p7|) for the 
kernel, it is possible to find the critical parameters that 
control the travelling waves for the mean- field limit {i.e. 
the infinite iVsat limit) of the evolution with saturation. 
A numerical solution of the condition jx'il) — xil) (fo'" 
A = 1 and B = 2) leads to 



7e = 0.69216--- , x(7c) = 5.2314- ■• 

and x'ilc) = 7.5582 • • • . 



(29) 



Expected properties of the amplitudes in these 
models 



At this stage, we can derive the properties of the trav- 
elling waves. The first step is to write down the BFKL 
evolution of the dipole densities. If Ui is the number of 
dipoles at lattice site i, we have 

^ (26) 

= ^ AlogBB(^"')'^nj -h^(B^ - l)nj. 

j>i j<i 

The amplitude T is just ni/Nsa±. 

The equation that corresponds to the BK equation, 
which preserves unitarity, is obtained from the latter by 
taking the minimum of n as given by the linear evolution 
and TVsat, at each rapidity step and for each dipole size. 

The FIP model has the same number of variables as 
the FKPP equation: one evolution variable and one spa- 
tial variable (which is the dipole size) in which diffusion 
can occur. It is a branching diffusion process in the dis- 
cretised space of dipole sizes, with a saturation condition 
that limits the number of objects (dipoles) of each size. 
From a general analysis of such models, we know that 
the realisations of the evolution of the scattering ampli- 
tude are noisy travelling waves whose properties may be 
obtained from the branching diffusion kernel. 

The eigenfunctions of the kernel of the fixed-impact 
parameter evolution are of the form B~'^'^ with the cor- 
responding eigenvalues 

X(7) = ^IWogBB-''^ B''^'^ +Y.^B'^ ^ 1) ^''^'^ 



k<0 



k<0 



AlogS 
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1 _^(7~1)A ^7A _ l' 



(27) 



Note that in the continuous limit obtained by letting A go 
to 0, one recovers the standard coUinear approximation 
to the BFKL kernel: 



7 1 -7 



(28) 



which confirms that we have kept the most important 
singularities. Furthermore, one may verify that the I/7 
and 1/(1 — 7) singularities are present independently of 



The travelling waves proceed to larger values of i as rapid- 
ity is increased. They are characterised by their position 
Ps, related to the saturation scale in QCD, and by the 
shape of their leading edge. Note that there is a certain 
freedom in the definition of the position of the front: we 
will specify our choice in Sec. IIV Al The large-rapidity 
velocity of the wave reads 



dps 
dY 



1 



logs 



x'ilc 



while the shape is given by 

T{i) ~ B-7c(A^-p,) 



(30) 



(31) 



in the forward part of the wave front defined by lA ^ ps. 
For a finite A'sat, the velocity is in general lower, and 
differs from the the mean-field velocity by amount of 
order 1/log iVgat- Furthermore, the front position ac- 
quires a dispersion from event to event, which goes like 
Y/ log'^ A^sat asymptotically for large Y and A^sat ■ 

If the property that all impact parameters evolve in- 
dependently of each other is true, travelling waves with 
the same properties are expected for T{i,b) in the AIP 
model, at each value of the impact parameter. This is 
what we shall investigate numerically in the last section. 



IV. NUMERICAL RESULTS 



We take as an initial condition a number A"sat of dipoles 
of size 1 (i = 0), uniformly distributed in impact pa- 
rameter between — ro/2 and ro/2 in the case of the AIP 
model. (This means that the impact parameters of the 
initial A'sat dipoles are randomly chosen according to a 
fiat distribution^). We vary A'sat between 10 and 200, 
and the cutoff k between 10~^ and 10~^. The impact pa- 
rameters bj that are considered are respectively 0, 10~^, 
10-4, 10-2 and 10-^ We set A = 1 and B = 2 in all 
that follows. 



^ We have checked that starting with a fixed initial distribution of 
dipoles in impact parameter leads to the same asymptotic results. 
See e.g. Fig. |H] further in this Section where the simulation for 
Naat = 100 up to y = 4 has been obtained by taking A'sat 
dipoles at 6 = as an initial condition: only differences in the 
pre-asymptotic regime are observed. 
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FIG. 2: Shape of the traveUing-wave front in the AIP model 
at central impact parameter, for three rapidities {Y — 1, 3 
and 5). Five different events are shown in dashed lines, while 
the average over many events is displayed with continuous 
lines. Nsat = 25 and k = 10~^. 



The number of events generated is typically 10** for 
each set of parameters, which allows one to measure the 
average and variance of the position of the travelling 
waves to a sufficient accuracy. 

The numerical data presented below were obtained us- 
ing two independent implementations of the Monte-Carlo 
event generator, which gave fully consistent results. 



A. Amplitudes at a given impact parameter 

For the needs of the discussion in this section, conve- 
nient values of A'sat and k are 25 and 10~^ respectively. 
The latter will be the default values, unless stated explic- 
itly. 

First, we can observe very clearly the propagation of 
the travelling waves, in the FIP model (not pictured) as 
well as in the AIP model at different impact parameters 
(see Fig. [S] similar curves would be obtained in the FIP 
model). Note that the different realisations shown in 
Fig. [2] sometimes overshoot 1. This is due to the fact 
that the unitarity constraint is only imposed at a finite 
number of points (in practice only 3, in order to save CPU 
time). Adding more points reduces the excess but does 
not change the final results for the wave front velocity 
and dispersion. 

To explore quantitatively the properties of the waves, 
we define their position ps as the largest i {i.e. the small- 



FIG. 3: Instantaneous velocity of the front in the FIP model, 
for A'sat ~ 25, 100, 1000 as a function of the rapidity. The 
mean-field value of the velocity x'ilc)/ log B is also shown in 
full line. 



est dipole size) for which T > 1/2, namely 

A^sat 



and iVsat X r(ps + l,fcj) < 



N,, 



(32) 



In this section, we will be interested in measuring the 
y-slope of the average of ps over events, as well as the 
variance of the front position 



a\Y) = (pUY)) - {ps{Y)y^ 



(33) 



First, let us address the FIP model in which there is 
no impact-parameter dependence. We plot the velocity 
of the travelling wave for different values of iVgat in Fig. [31 
and the corresponding variance in Fig. |4l both as a func- 
tion of the evolution variable Y. We see that the velocity 
reaches quite quickly its asymptotic regime: it becomes 
stable as soon as y > 2, for any value of j. The same 
holds for the growth rate of the variance. A fit in the 
range 2 < Y < 5 can safely be extrapolated to y > 5. 
(This remark will be useful in the case of the AIP model, 
where we have to restrict the calculation to F < 5 for 
technical reasons.) The results of the fits are shown in 
Tab. [H for several values of iVgat- 

Now we repeat the same calculation in the AIP model, 
for y < 5. The choice of the upper limit on Y is dic- 
tated by numerical limitations: if Y gets larger, then the 
travelling waves are likely to explore too small dipole sizes 
which cannot be resolved numerically for lack of accuracy 
(see Sec. IIII A)) . We first check that the front velocity is 



{pl(Y)) - (psiY))-' 



{ps{Y + dY,b)-ps(Y,b})/dY 





FIG. 4: Variance of the front position as a function of the 
rapidity in the FIP modeL Three different values of A'sat are 
displayed. A linear fit is indistinguishable from the numerical 
data. 



FIG. 5: Instantaneous velocity of the front in the AIP model, 
for A'sat ~ 25 and at different impact parameters. The curves 
are all almost perfectly superimposed. 
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10* 10360 28365 
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25 


7.33 


8.26 


8.38 


4.62 


3.07 


3.06 
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50 


7.84 


8.60 
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3.76 


2.44 


2.50 


10* 35127 7102 


100 


8.24 


8.87 
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3.11 


1.99 


2.25 


10* 15851 3735 


100* 


— 


— 


9.03* 


— 


— 


2.03* 


— — 10** 


200 


8.57 


9.11 


— 


2.57 


1.73 


— 
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TABLE I: Velocity and rate of growth of the dispersion of 
travelling waves in the FIP and AIP model, for different values 
of A'sat and k (in the case of the AIP model). For Asat = 100 
and K = 10-"^, we have also performed a calculation up to 
y = 4 only (data denoted by a star) in order to be able to 
collect enough statistics. 



independent of the impact parameter that is considered, 
Fig. [5] As this is indeed the case, we are free to focus 
on central impact parameters. Figure [S] shows travelling 
waves averaged over many events at central impact pa- 
rameter, for different values of the cutoff k. (Note that 
the shape of the front in its forward part is consistent 
with the theoretical expectation in Eq. ([7])). Figures [7] 
and [H] show the average of the instantaneous velocity and 
of the dispersion of the front position as a function of the 
rapidity. Again, a stationary velocity is reached after 1 
to 2 units of rapidity. As soon as this regime is reached, 
the dispersion starts to scale linearly in Y , as expected. 



The diffusion coefficient D = da'^ /dY, which is the slope 
of this line, quantifies the amount of fluctuations. The 
values of the velocity and of the diffusion coefficient are 
reported in Tab. HI also for different values of A'sat (and 
k) . The numbers have been extracted from a fit between 
Y — 2 and 5 for the velocity, and between Y = 2 and 
4.5 for the diffusion coefficient. We have limited the fit 
range for the latter* to avoid edge effects related to the 
fact that the far tail of the front sometimes reaches our 
cut-off on dipole sizes (see e.g. Fig. [S]), leading to the 
small turnover for the dispersion observed in Fig. [51 

From Figs. [7] and [H we already see that there is no 
strong dependence on the value of k. To further study 
this dependence, we have computed the asymptotic ve- 
locity for different values of k and A'sat = 25. One sees 
in Fig. [9] a saturation for k < 10"'^ and a very small de- 
pendence on K up to ~ 10^^, confirming that the cutoff 
K does not strongly affect the observables. 

There is an important difference between AIP and FIP, 
where similar results were expected if the hypothesis of 
local evolution in impact parameter were literally correct. 
Instead, comparing Figs. [3] and [5] with Figs. [7] and [51 we 
observe that the FIP model has more fluctuations than 
the AIP model, as if the effective number of particles were 
lower, or if some averaging were effectively performed in 
the AIP evolution. In order to enforce the matching of 
the velocities and diffusion coefficients of the two models. 



* Since the velocity instead depends on what happens close to the 
saturation scale, the fit can safely be extended to Y = 5. 
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FIG. 6: Average of the event-by-event amplitude in the AIP 
model as a function of the dipole size at zero impact param- 
eter, for three different rapidities. Two values of the cutoff k 
are compared and A^sat is set to 25. 



FIG. 7: Instantaneous velocity of the front in the AIP model, 
for A'^sat ~ 25 and A'sat ~ 100, and different values of the 
cutoff K. The impact parameter b is set to 0. The mean-field 
result is also represented. 



we should take a value of Nsat of about 100 for the FIP 
model, i.e. about 4 times that of the AIP model (with 
iVsat = 25). We will try and investigate more systemat- 
ically this discrepancy below. First, let us examine the 
correlations between travelling waves at different impact 
parameters. 



invert the above formula for the relevant values of 62~&l)• 
These rapidities are denoted by a cross in Fig. [TO] for the 
considered impact parameter differences. Our numeri- 
cal results for the correlations are nicely consistent with 
this estimate, since the correlations start to saturate to a 
constant value precisely on the right of each such cross. 



B. Correlations in impact parameter 

Fig. [TO] represents the correlations between the posi- 
tions of the wave fronts at different impact parameters 
in the AIP model, defined as 



{p,{Y,h)pAY,b2)) - {ps{Y,bi)){p,(Y,b2)). 



(34) 



Once again, the different choices of k lead to very similar 
results. 

We also see very clearly the successive decouplings of 
the different impact parameters, from the most distant to 
the closest one, as rapidity increases. Indeed, the corre- 
lation functions flatten after some given rapidity depend- 
ing on the difference in the probed impact parameters, 
which means that the evolutions decouple. This decou- 
pling is expected as soon as the travelling wave front 
reaches dipole sizes which are smaller than the distance 
between the probed impact parameters, i.e. at Y such 
that I62 - bi\ « 1/QsiY) = B-P-^^\ From the data 
for ps{Y)^ we can estimate quantitatively the values of 
the rapidities at which the travelling waves decouple be- 
tween the different impact parameters. (It is enough to 



A'^sat -dependence and difference betw^een AIP 
and FIP 



We have observed that the travelling wave fronts in the 
AIP model at a given impact parameter look very much 
like those appearing in the FIP model, but the parame- 
ter A'sat should be changed by some factor to allow for 
a matching of their characteristics. In this section, we 
investigate more quantitatively this point. 

To this aim, we compute the front velocity V and diffu- 
sion coefficient D = da^/dY in the AIP model at central 
impact parameter as a function of iVgat ■ We use different 
values of the cutoff k (in practice k — 10~^ and 10~^). A 
similar calculation can be performed for the FIP model 
(in this case, there is no cutoff k to be considered). The 
results are displayed in Fig. [TlJ 

We notice that the shapes of these curves look similar, 
in the FIP and AIP model for V as well as for D, except 
maybe for very small values of Ngat (10 and 15) . However, 
we may superpose the different curves only at the price 
of performing a rescaling of A'sat (which corresponds to a 
shift of the curves on the logarithmic scale chosen in the 
plot of Fig.fTTjl. Empirically, we find that the AIP model 
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FIG. 8: The same for the variance of the position of the front 
at central impact parameter. Note that for A'sat = 100 and 
K = 10"'^, we have data only up to Y = 4. 
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FIG. 9: K dependence of the asymptotic front velocity for 
iVsat = 25. 



is equivalent to the FIP model at each impact parameter 
if one sets A'sat for the latter about 3.8 times larger than 
for the former. 

We may propose an interpretation of the observed dis- 
crepancy. In the course of the evolution, travelhng waves 
at different impact parameters communicate by exchang- 
ing dipoles because dipole splitting is not strictly local in 
coordinate space: starting from any given step, the chain 
of subsequent splittings extends over an area that is sig- 



FIG. 10: Correlations between the front position at different 
impact parameters as a function of the rapidity for A'sat = 25. 
The crosses denote the rapidity Y at which the (inverse of the) 
saturation scale iJ~''=(^) coincides with the distance in impact 
parameter space between the probed points. 



nificantly larger than the area set by the initial dipole. 
In this way, the waves may equalise their velocities, and 
this leads to an extra dynamical averaging with respect 
to FIP models, at any impact parameter, continuously 
at any step in the evolution. This effective averaging re- 
sults in a reduction of the fluctuations in the positions 
of the waves, and thus to an increase of the front ve- 
locity, which becomes closer to the mean-field velocity, 
and in particular different from the expected velocities in 
the corresponding one-dimensional travelling wave model 
(FIP model). 

As we have observed in Fig. [Til the scaling factor be- 
tween A'sat and the effective maximum number of dipoles 
in the AIP model seems to be independent of A'sat ■ How- 
ever, we also see that a factor log(A'sat)/7c (the length of 
the front) between AIP and FIP is also consistent with 
our numerical observation. Unfortunately, we have not 
been able to find a satisfactory analytic explanation for 
this factor. 

In any case, since we find that the discrepancy amounts 
to a simple rescaling of A'sat by a constant or a slowly 
varying factor (at least in the range of the values of A'sat 
that we were able to explore) , it is likely that it is a sub- 
leading effect, and that the AIP and FIP models would 
agree asymptotically for large A'sat- Hence we do not 
think that this slight mismatch between the AIP and 
the FIP models spoils in any way the one-dimensional 
reaction-diffusion picture of high-energy scattering. 
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FIG. 11; Steady-regime velocities (left) and diffusion coefficients (right) of the wave front in the FIP model (full lines) and at 
central impact parameter in the AIP model (dashed and dotted lines). The two leftmost curves in each plot correspond to the 
actual results of the calculation; The two rightmost ones are the same but shifted. The shift is a scaling factor of A^sat ~ I/q^ 
equal to 3.8. The curves with filled black squared correspond to the AIP model shifted by — logAgat- For Agat ~ 100 and 



10" 



we have performed a fit up to F = 4 only (see Tab. |T]|. 



V. CONCLUSION 

One condition for high-energy evolution in QCD to be 
in the universality class of a one-dimensional reaction- 
diffusion processes is that the partonic evolution remains 
local in impact-parameter space. In this paper, with the 
help of a toy model (called "AIP" ) engineered to incor- 
porate the essential features of QCD and to be man- 
ageable numerically, we have checked that this is indeed 
the case. We have measured the properties of travelling 
waves at given impact parameters, and showed that the 
correlations between waves at different points in coordi- 
nate space vanish when the rapidity gets large enough, 
as expected for a local evolution. 

Furthermore, we have compared the properties of the 
travelling waves that appear in this model with the fea- 
tures of their counterparts in a model that has no impact 
parameter dependence (called "FIP"), but which is iden- 
tical to the former in all other respects. We have found 
that the results are similar, except that the number of 
partons at saturation should be larger by a factor of the 
order of 4 in the FIP model, suggesting that the latter has 
more fluctuations. It is also plausible that this factor has 
a slow A'sat-dependence and may vary like — log A^sat- At 
this stage, our numerical data do not enable us to check 
whether one of these guesses is literally correct. 

We do not have a final explanation for the discrep- 
ancy between the AIP and FIP models yet. However, 
the fact that the factor by which A'sat should be scaled 
for the matching to occur does not exhibit a strong A^sat- 



dependence (remember that Agat ~ l/o^s ii^ QCD) leads 
us to think that it is a subleading effect in the large- 
A^sat (small-as) limit, and thus, that the one-dimensional 
reaction-diffusion picture is indeed valid, up to the rele- 
vant replacement of the parameter A^sat by some effective 
number of particles relevant to the propagation of the 
travelling waves, yet to be understood more precisely. 



Note that, though a full study with two transverse de- 
grees of freedom would be of great interest, we believe 
that our one-dimensional picture grasps the important 
aspects of the problem and, based on universal properties 
of the reaction-diffusion systems, we expect our results to 
hold for fuU QCD. 



Finally, we can think of several other improvements 
beyond our study. We have considered only one satura- 
tion mechanism, which consisted in merely vetoing par- 
ton splittings in phase space cells already occupied by 
at least A"sat particles. We could try and check whether 
other saturation mechanisms (such as parton recombi- 
nation) would bring about more correlations in impact 
parameter. It would also be desirable to achieve some 
analytical understanding, in particular of the rescaling 
factor between the AIP and FIP model. To this aim, it 
may be interesting to study the chains of successive dipole 
splittings in the AIP model, in order to find out whether 
rare fiuctuations spread in a wider-than-expected area in 
impact parameter. 
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